home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / vel.pro < prev    next >
Text File  |  1997-07-08  |  5KB  |  176 lines
  1. ; $Id: vel.pro,v 1.4 1997/01/15 03:11:50 ali Exp $
  2.  
  3. function vel_mybi,a,x,y
  4. on_error,2                      ;Return to caller if an error occurs
  5. sizea=size(a)
  6. nx=sizea[1]
  7. i=long(x)+nx*long(y)
  8. q=y-long(y)
  9. p=x-long(x)
  10. q1 = 1.-q
  11. p1 = 1.-p
  12.  
  13. ; Weighting factors were wrong for a(i+1) & a(i+nx), switched them.
  14.  
  15. aint=p1*q1*a[i] + p*q1*a[i+1] + q*p1*a[i+nx] + p*q*a[i+nx+1]
  16. return,aint
  17. end
  18.  
  19. PRO ARRHEAD,X
  20.  
  21. ON_ERROR,2                      ;Return to caller if an error occurs
  22. theta = 30 * !radeg
  23. TANT = TAN(THETA)
  24. NP=3.0
  25. SCAL=8.
  26.  
  27. SX=SIZE(X)
  28. N=SX[2]
  29.  
  30.  
  31. BIGL=SQRT((X[*,N-4,0]-X[*,N-5,0])^2+(X[*,N-4,1]-X[*,N-5,1])^2)
  32. wbigl=where(BIGL ne 0.0)
  33. wnbigl=where(bigl eq 0.0, count)
  34. LL  = SCAL*TANT*BIGL[wbigl]/NP
  35.  
  36. DX = LL*(X[wbigl,N-4,1]-X[wbigl,N-5,1])/BIGL[wbigl]
  37. DY = LL*(X[wbigl,N-4,0]-X[wbigl,N-5,0])/BIGL[wbigl]
  38.  
  39. XM = X[wbigl,N-4,0]-(SCAL-1)*(X[wbigl,N-4,0]-X[wbigl,N-5,0])/NP
  40. YM = X[wbigl,N-4,1]-(SCAL-1)*(X[wbigl,N-4,1]-X[wbigl,N-5,1])/NP
  41.  
  42. X[wbigl,N-3,0] = XM-DX
  43. X[wbigl,N-2,0] = X[wbigl,N-4,0]
  44. X[wbigl,N-1,0] = XM+DX
  45.  
  46. X[wbigl,N-3,1] = YM+DY
  47. X[wbigl,N-2,1] = X[wbigl,N-4,1]
  48. X[wbigl,N-1,1] = YM-DY
  49.  
  50. if count ge 1 then begin  ;No head for 0 length
  51.     X[wnbigl,N-3,0] = x[wnbigl,n-4,0]
  52.     X[wnbigl,N-2,0] = X[wnbigl,n-4,0]
  53.     X[wnbigl,N-1,0] = X[wnbigl,n-4,0]
  54.     
  55.     X[wnbigl,N-3,1] = X[wnbigl,N-4,1]
  56.     X[wnbigl,N-2,1] = X[wnbigl,N-4,1]
  57.     X[wnbigl,N-1,1] = X[wnbigl,N-4,1]
  58.     endif
  59.  
  60. return
  61. END
  62.  
  63. function arrows,u,v,n,length,nsteps=nsteps
  64. on_error,2                      ;Return to caller if an error occurs
  65. su=size(u)
  66. nx=su[1]
  67. ny=su[2]
  68.  
  69. lmax=sqrt(max(u^2+v^2))        ;Max vector length
  70. lth=1.*length/lmax/nsteps
  71. xt=randomu(seed,n)        ;Starting position
  72. yt=randomu(seed,n)
  73. x=fltarr(n,nsteps+3,2)
  74. x[0,0,0]=xt
  75. x[0,0,1]=yt
  76. for i=1,nsteps-1 do begin
  77.  xt[0]=(nx-1)*x[*,i-1,0]
  78.  yt[0]=(ny-1)*x[*,i-1,1]
  79.  ut=vel_mybi(u,xt,yt)
  80.  vt=vel_mybi(v,xt,yt)
  81.  x[0,i,0]=x[*,i-1,0]+ut*lth
  82.  x[0,i,1]=x[*,i-1,1]+vt*lth
  83. end
  84. ARRHEAD,X
  85. return,x<1.0>0.0
  86. end
  87.  
  88.  
  89.  
  90. PRO VEL,U,W,LENGTH=length,XMAX=xmax, nvecs = nvecs, nsteps = nsteps, $
  91.     title = title
  92. ;+
  93. ; NAME:
  94. ;    VEL
  95. ;
  96. ; PURPOSE:
  97. ;    Draw a velocity (flow) field with arrows following the field 
  98. ;    proportional in length to the field strength.  Arrows are composed 
  99. ;    of a number of small segments that follow the streamlines.
  100. ;
  101. ; CATEGORY:
  102. ;    Graphics, two-dimensional.
  103. ;
  104. ; CALLING SEQUENCE:
  105. ;    VEL, U, V
  106. ;
  107. ; INPUTS:
  108. ;    U:    The X component at each point of the vector field.  This 
  109. ;        parameter must be a 2D array.
  110. ;
  111. ;    V:    The Y component at each point of the vector field.  This 
  112. ;        parameter must have the same dimensions as U.
  113. ;
  114. ; KEYWORD PARAMETERS:
  115. ;    NVECS:    The number of vectors (arrows) to draw.  If this keyword is
  116. ;        omitted, 200 vectors are drawn.
  117. ;
  118. ;    XMAX:    X axis size as a fraction of Y axis size.  The default is 1.0.
  119. ;        This argument is ignored when !p.multi is set.
  120. ;
  121. ;    LENGTH:    The length of each arrow line segment expressed as a fraction 
  122. ;        of the longest vector divided by the number of steps.  The 
  123. ;        default is 0.1.
  124. ;
  125. ;    NSTEPS:    The number of shoots or line segments for each arrow.  The
  126. ;        default is 10.
  127. ;
  128. ;    TITLE:    A string containing the title for the plot.
  129. ;    
  130. ; OUTPUTS:
  131. ;    No explicit outputs.  A velocity field graph is drawn on the current
  132. ;    graphics device.
  133. ;
  134. ; COMMON BLOCKS:
  135. ;    None.
  136. ;
  137. ; SIDE EFFECTS:
  138. ;    A plot is drawn on the current graphics device.
  139. ;
  140. ; RESTRICTIONS:
  141. ;    none
  142. ;
  143. ; PROCEDURE:
  144. ;    NVECS random points within the (u,v) arrays are selected.
  145. ;    For each "shot" the field (as bilinearly interpolated) at each
  146. ;    point is followed using a vector of LENGTH length, tracing
  147. ;    a line with NSTEPS segments.  An arrow head is drawn at the end.
  148. ;
  149. ; MODIFICATION HISTORY:
  150. ;    Neal Hurlburt, April, 1988.
  151. ;    12/2/92    - modified to handle !p.multi (jiy-RSI)
  152. ;       7/12/94 HJM - Fixed error in weighting factors in function
  153. ;                     vel_mybi() which produced incorrect velocity vectors.
  154. ;
  155. ;-
  156.  
  157. on_error,2                      ;Return to caller if an error occurs
  158. if n_elements(Nvecs) le 0 then nvecs=200
  159. if n_elements(nsteps) le 0 then nsteps = 10
  160. if n_elements(length) le 0 then length=.1
  161. if n_elements(title) le 0 then title='Velocity Field'
  162. X=ARROWS(U,W,Nvecs,LENGTH, nsteps = nsteps)
  163.  
  164. if (!p.multi[1] eq 0 and !p.multi[1] eq 0) then begin
  165.    if (n_elements(xmax) eq 0) then xmax = 1.0
  166.    IF XMAX GT 1. THEN position=[0.20,(0.5-0.30/XMAX),0.90,(0.5+0.40/XMAX)]$
  167.       else position=[(0.5-0.30*XMAX),0.20,(0.5+0.40*XMAX),0.90]
  168.    plot,[0,1,1,0,0],[0,0,1,1,0],title=title,pos=position
  169. endif else begin
  170.    plot,[0,1,1,0,0],[0,0,1,1,0],title=title
  171. endelse
  172.  
  173. FOR I=0,Nvecs-1 DO PLOTS,X[I,*,0],X[I,*,1]
  174. RETURN
  175. end
  176.